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I .  INTRODUCTION 


Liquid  regenerative  propellant  guns  have  been  proposed  as  a 
possible  substitute  for  the  traditional  solid  propellant  gun.  Morrison 
et.  al.  have  covered  the  background  and  possible  applications  of  liquid 
propellant  guns . *  Figure  1  shows  the  design  for  a  typical  regenerative 
gun.  The  monopropellant  is  originally  in  the  liquid  reservoir.  A 
primer  pressurizes  the  combustion  chamber.  This  starts  the  regenerative 
piston  moving.  Liquid  is  injected  between  the  moving  piston  and  the 
fixed  central  bolt.  As  the  liquid  ignites  and  burns,  the  rise  in 
pressure  accelerates  the  piston  motion  and  the  subsequent  fluid 
injection.  The  taper  at  the  back  of  the  central  bolt  slows  down  the 
piston  at  the  end  of  the  injection  stroke. 


Existing  liquid  propellant  gun  codes^"5  involve  a  number  of 
simplifying  assumptions.  In  this  paper  the  modeling  of  the  injection  of 
the  liquid  through  the  orifice  between  the  piston  and  the  bolt  is 
described.  In  present  gun  codes,  this  is  usually  approximated  by 
steady- state  Bernoulli  flow.  The  discharge  coefficient  is  constant 
throughout  the  firing  cycle.  This  constant  value  is  often  chosen  so  as 
to  match  a  desired  maximum  pressure.  However,  analysis  of  data  from  a 
30 -mm  regenerative  gun  fired  at  the  BRL  indicate  that  a  steady  discharge 
coefficient  is  not  a  good  approximation.^"'7  When  cast  in  terms  of 
steady- state  Bernoulli  flow,  the  discharge  coefficient  starts  small  and 
only  gradually  builds  up  to  something  close  to  a  steady  state  value. 

The  transient  lumped  parameter  model  (LPIN)  was  described  in  the  latter 


paper,  but  the  model  did  not  agree  well  with  the  experimental  data. 


Science  Applications  International  Corporation  (SAIC)  has  a 
contract  to  develop  a  two-dimensional  model  for  the  injection  process. 
The  model  is  based  on  an  existing  code,  MAGIC  (Modeling  Algorithm  of  a 
Generic  Internal  Combustion  Engine),  developed  to  model  injection  and 


combustion  in  a  diesel  engine.  This  has  been  modified  to  model  the 


liquid  reservoir  in  a  regenerative  liquid  propellant  gun.  In  addition, 
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the  model  has  been  used  in  the  analysis  of  interior  ballistic  processes 
of  bulk  loaded  liquid  propellant  guns . ^ 

11  1? 

The  MAGIC  code  is  based  on  algorithms  developed  at  Los  Alamos. 

It  is  a  two-dimensional  time-marching  code  that  solves  finite-difference 
approximations  of  the  partial  differential  equations  describing 
compressible  fluid  dynamics.  A  control  volume  approach  has  been  used 
with  a  staggered  grid.  The  grid  can  vary  quite  generally  in  time  (ALE  - 
arbitrary  Lagrangian-Eulerian  mesh) .  Two  temporal  differencing  options 
are  provided.  The  first  option  is  an  explicit  difference  scheme  whose 
primary  stability  limitation  is  the  Courant- Friedrichs -Levy  condition. 

A  partially  implicit  scheme  can  also  be  used  (ICE  method).  Centered 
spatial  differencing  is  utilized  for  all  stress  terms  in  the  equations. 
The  spatial  advective  terms  can  be  approximated  by  purely  centered  to 
full  donor  cell  differencing. 

At  present  there  are  problems  with  this  code.  The  primary 
difficulty  is  that  the  code  is  based  on  a  first  order  finite  difference 
approximation.  The  convergence  cf  such  a  code  is  very  slow.  Since 
discharge  coefficients  must  be  predicted,  very  accurate  solutions  are 
required.  Even  on  the  BRL  CRAY-XMP,  a  sufficiently  fine  grid  cannot  be 
run,  because  of  both  memory  and  time  restrictions . 

A  secondary  problem  is  that  the  adaptive  gridding  algorithm  is  not 
general  enough.  The  MAGIC  code  was  designed  for  use  in  a  simple 
cylinder  with  one  end  moving.  The  motion  of  a  complicated  piston  shape 
moving  past  a  tapered  bolt  is  more  difficult.  Currently,  SAIC  is 
modifying  the  code  to  overcome  these  problems. 

Because  of  these  complexities  a  one -dimensional  injection  code 
(ODIN)  has  been  developed.  The  algorithm  is  simple  enough  that  the 
complicated  piston  and  bolt  shapes  can  be  readily  implemented,  and 
numerical  convergence  can  be  easily  obtained.  Even  after  the  two- 
dimensional  code  is  fully  developed,  the  one -dimensional  code  will  be 
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useful  for  making  preliminary  runs . 

This  report  discusses  the  development  of  the  one  dimensional  code. 
For  a  simple  test  problem,  the  code  is  compared  with  the  MAGIC  code. 

The  results  are  quite  similar.  When  applied  to  the  30-mm  gun  problem, 
the  one -dimensional  code  still  gives  poor  agreement  with  the  experiment. 
It  does  agree  closely  with  an  earlier  lumped  parameter  injection  model 
(LPIN) . 

The  results  suggest  that  even  the  two-dimensional  model  will  not  be 
able  to  predict  the  injection  process.  There  is  apparently  some 
phenomenon  we  are  not  properly  taking  into  account. 

II.  GOVERNING  EQUATIONS 

The  fundamental  governing  equations  are  conservation  of  mass  and 
momentum  in  the  liquid  propellant.  Preliminary  runs  with  the  MAGIC  code 
indicate  that  the  energy  increase  in  the  liquid  due  to  the  work  of 
compression  is  negligible,  so  the  energy  equation  is  not  implemented.  A 
control  volume  approach  is  used.  The  equations  are  expressed  in 
integral  from  for  a  control  element  which  may  be  in  motion  with  an 
arbitrary  prescribed  velocity  U. 

The  continuity  equation  is^ 

3/at  J  pdv  -  I  p(U  -  u) •  n  dS  -  0  (1) 

V(t)  JS(t) 

where  V  is  the  volume  of  the  element,  S  is  the  surface  area,  &  is  the 
outward  normal  vector,  p  is  the  density,  and  u  is  the  fluid  velocity. 
Although  the  liquid  is  almost  incompressible,  over  the  large  pressure 
ranges  of  interest  the  density  changes  will  be  important.  The  momentum 
equation  is^® 


I 


’’WjK 


v.v.-.v-  /,  . 


d/dt  J  pudV  -  I  pu(U  -  u)  •  ndS  +  J  pndS 

V(t)  S(t)  S(t) 

-  J  u  Vu  ndS  -  0  (2) 

S(t) 

where  p  is  the  pressure  and  p  is  the  dynamic  viscosity.  The  kinematic 
viscosity  has  been  measured  for  several  propellants  as  a  function  of 
temperature.  In  this  article,  HAN1846  at  room  temperature  will  be 
modeled  for  which  the  kinematic  viscosity  is  0.05  stokes.  The  dynamic 
viscosity  p  —  pi/ . 

The  equation  of  state  is^'^ 

p  -  K1/K2  t (P/P0  )K2  -  1]  (3) 

where  is  the  bulk  modulus  at  zero  pressure  and  K2  is  the  derivative 
of  the  bulk  modulus  with  respect  to  pressure  (bulk  modulus  K  -  +  K2 

p).  For  HAN1846 ,  pQ-1.43  g/cc,14  K:  -  5350  MPa,  and  K2  -  9. 11. 15 

For  the  initial  conditions  the  fluid  velocity  is  chosen  to  be 
identically  zero  and  the  pressure  field  to  be  uniform.  Other  initial 
conditions  have  been  tried,  but  if  they  are  not  physically  reasonable, 
large  oscillations  are  set  up  in  the  solution. 

The  boundary  conditions  will  be  chosen  for  physical  reasons.  In 
the  problem  of  interest,  we  know  the  piston  velocity  and  combustion 
chamber  pressure.  As  seen  below,  this  will  translate  into  specifying 
the  influx  velocity  and  the  outflux  pressure. 

III.  FINITE  DIFFERENCE  APPROXIMATION 


The  finite  difference  scheme  is  taken  from  the  MAGIC  code.  A 
staggered  grid  is  used  (see  Figure  2).  Vector  quantities  are  defined  on 
the  main  grid  (solid  lines) .  Scalar  quantities  are  defined  at  points 
midway  between  the  main  grid  values  (dotted  lines).  For  flexibility  in 
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implementing  boundary  conditions,  the  grid  is  extended  past  the  physical 
boundary.  Let  N-l  be  the  number  of  intervals  in  the  region  of 
integration.  The  vector  grid  then  goes  from  J-l  to  N+l.  J-l  is  the 
left  boundary,  and  J-N  is  the  right  boundary.  The  scalar  grid  goes  from 
1-1  to  N.  The  grid  point  I  is  midway  between  the  vector  points  J  and 
J+l.  Unless  otherwise  stated,  it  is  assumed  that  I-J.  Two  different 
letters  are  used  to  distinguish  between  the  vector  grid  and  the  scalar 
grid. 

Now  consider  Eq.  (1),  the  continuity  equation,  at  point  I.  The 
control  volume  for  scalar  quantities  is  taken  to  be  from  point  J  to  J+l. 
Let  Mj  be  the  mass  in  the  control  volume  at  the  present  time,  and  Mjn  be 
the  mass  at  the  next  time  step.  Note  that  tyj  -  Vj .  The  change  in 
mass  is  due  to  mass  flux  into  the  control  volume  minus  mass  flux  out. 
Using  forward  time  differences, 


(M?  '  Mi)/At  “  Pj<uj  *  UJ>AJ  -  *\j+l(uj+l  *  uJ+l)AJ+l 

where  A  is  the  cross  sectional  area.  The  density  is  not  known  at  the 
vector  grid  points,  so  it  must  be  approximated.  To  improve  the 
stability  of  the  scheme,  upwind  differencing  is  used.^  The  physical 
rationale  is  that  the  density  at  a  point  will  be  primarily  determined  by 
the  upstream  conditions.  So  let 

i Pi . 1  u j  >  0 

Pj  ‘  n  (5) 

l  pi  uj  <  0 

The  density  at  the  new  time  step  is  obtained  by  dividing  the  mass  by  the 
control  volume . 

Next  consider  Equation  (2),  the  momentum  equation,  at  point  J.  The 
control  volume  for  the  vector  quantities  is  taken  from  point  1-1  to  I. 
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KMjuj)n  -  mjuj]/^  -  -80(Pl  -  Pl-1>(AI  +  AI-1>/2 

+  »/[pjAjduj  -  +  ^j(Aj  -  Aj_^)duj] 

+  ^I-1AI-1UI-1^UI-1  *  ^i-i)  "  PjAjUj(uj  -  Uj)  (6) 

The  pressure  term  includes  not  only  the  pressure  at  the  left  and  right 
sides  of  the  control  volume,  but  also  the  pressure  on  the  side  walls 
(area  -  Aj  -  Aj_^).  This  of  course  only  has  an  effect  when  the  area  is 
changing.  The  pressure  on  the  side  walls  is  taken  to  be  the  average  of 
the  pressures  at  the  boundaries  of  the  control  volume.  The  resulting 
three  components  are  combined  into  the  above  form.  Similarly,  the 
viscosity  term  has  three  components.  The  velocity  derivatives  are 
approximated  by  central  differences . 


IS 


duI  -  (uj  -  uJ.1)/(xJ  -  xj.p  (7) 

duj  “  <uJ+l  '  uJ-l)/<xJ+l  '  XJ-1> 

The  velocity  at  the  side  walls  is  taken  to  be  zero  (no  slip  condition) . 
The  convection  terms  consist  of  the  momentum  flux  into  the  control 
volume  minus  the  momentum  flux  out.  As  before,  upwind  differencing  is 
used  for  the  convection  term. 
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(9) 


The  mass  Mj  of  the  vector  control  volume  must  also  be  computed.  This  is 
a  weighted  average  of  the  masses  of  the  two  scalar  control  volumes. 


To  maintain  stability,  the  time  step  will  still  be  constrained  by 
the  Courant- Friedrich -Levy  (CFL)  restriction^ 
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where  a  is  the  speed  of  sound  in  the  liquid,  given  by 

a  “  J  g 0K/P 


(10) 


(11) 


IV.  IMPLEMENTATION 

The  program  was  developed  to  model  flow  through  an  annular  orifice . 
The  initial  shape  of  the  orifice  is  determined  by  the  central  bolt 
radius  and  the  piston  radius.  The  central  bolt  consists  of  straight 
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sections.  The  piston  contains  two  sections  that  are  arcs  of  circles. 

The  coordinate  system  is  attached  to  the  piston.  The  velocity  of  the 
back  wall  is  considered  known  (at  J-l)  .  The  back  wall  and  center  bolt 
move  past  the  piston.  In  the  actual  fixture,  the  piston  is  the  moving 
part.  For  a  one -dimensional  code,  it  is  irrelevant  which  part  is 
considered  fixed  and  which  part  is  considered  moving. 

To  specify  the  initial  grid,  the  number  of  intervals  (N-l)  and  an 
expansion  factor  FAC  are  entered.  A  grid  is  generated  such  that  the 
largest  interval  (at  the  left)  is  FAC  times  as  large  as  the  smallest 
interval  (at  the  right) .  The  grid  size  changes  linearly  between  the 
boundaries.  As  time  goes  on,  the  back  wall  moves  forward.  The  velocity 
of  the  first  grid  point,  at  the  back  wall,  is  known.  The  velocity  of 
the  right  hand  boundary  (J-N)  is  set  to  zero.  The  velocities  of  the 
other  grid  points  are  found  by  linear  interpolation.  The  intervals  at 
the  left  are  compressed  the  most,  which  is  why  they  are  chosen  to  be 
initially  larger.  As  the  grid  moves,  the  radii  and  cross  sectional 
areas  must  be  recomputed.  Normally,  a  value  of  FAC-2  was  chosen  for 
these  calculations. 


The  influx  velocity  (J-l)  and  the  outflux  pressure  (I-N)  are  user 
supplied  boundary  conditions.  The  pressure  and  density  at  the  influx 
boundary  (J-l)  and  the  velocity  at  the  outflux  boundary  (J-N+l)  are 
determined  simply  by  first-order  extrapolation.  This  has  been  shown  to 
be  as  stable  and  accurate  as  more  complicated  schemes  for  at  least  some 


problems . 
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The  time  step  is  based  on  the  CFL  condition,  equation  (7).  The 
Courant  number  C  is  computed  based  on  the  vector  control  volumes.  The 
velocity  u  is  taken  to  be  the  average  of  the  velocities  at  the 
boundaries  of  the  control  volume.  The  smallest  Courant  number  is 
multiplied  by  a  safety  factor  (usually  0.75)  to  obtain  the  next  time 
interval.  The  partially  implicit  pressure  iteration,  which  would  allow 
a  larger  time  step,  has  not  been  implemented  in  the  one -dimensional  code 
(see  discussion  below). 

The  time  step  integration  proceeds  in  two  main  stages.  First  a 
Lagrangian  calculation  is  performed,  assuming  that  the  grid  moves  with 
the  fluid.  The  convection  terms  do  not  contribute  to  this  part  of  the 
calculation.  Then  the  rezone  phase  takes  into  account  the  movement  of 
the  fluid  with  respect  to  the  grid  (convection  terms).  In  detail: 

1.  Compute  the  size  of  the  next  time  step  At  based  on  the  CFL 
condition. 

2.  Compute  the  new  velocity  for  the  back  wall  and  the  outflux  pressure 
(considered  known) . 

3.  Interpolate  to  get  the  velocity  of  all  the  grid  points. 

4.  Find  the  position  of  the  new  grid.  Assume  the  grid  velocity  over 
the  time  interval  is  the  average  of  the  old  and  the  new  velocity. 


5.  Compute  the  cross  sectional  areas  for  the  new  grid.  Compute  the 
control  volumes  for  the  new  grid. 


6.  Compute  the  new  mass  times  fluid  velocity  based  only  on  the  stress 
terms  (pressure  and  viscosity)  in  equation  (6).  Assume  for  the  moment 
that  the  grid  moves  with  the  fluid  (Lagrangian  calculation)  so  the 
convection  terms  do  not  contribute. 

7.  Compute  the  mass  in  the  vector  control  volume  by  a  weighted  average 
of  the  old  masses  in  the  scalar  control  volumes. 

8.  Compute  the  intermediate  fluid  velocities. 

9.  Compute  the  new  masses,  equation  (4),  using  the  updated  velocities. 
Find  the  new  densities  using  the  new  values  for  the  control  volumes. 

10.  From  the  equation  of  state  (3),  compute  the  new  pressures. 

11.  Using  the  new  densities  and  pressures,  update  the  mass  times 
velocity  based  on  the  convection  term  in  equation  (6). 

12.  Compute  the  mass  in  the  vector  control  volume  by  a  weighted  average 
of  the  new  masses  in  the  scalar  control  volumes. 

13.  Compute  the  final  fluid  velocities. 

The  results  can  be  summarized  in  terns  of  the  discharge 
coefficient.  The  steady  state  Bernoulli  equation  is 

Un  -  CD  J  2g0  (Pl  -  p n)/p  +  uL2  (12) 

where  p  is  the  average  density  in  the  orifice.  The  subscripts  1  and  N 
denote  the  boundaries.  This  equation  is  derived  from  the  steady  state 
momentum  equation  (ignoring  friction  losses).  The  discharge  coefficient 


Cjj  is  an  empirical  number  to  take  into  account  the  loss  terms  and  the 
error  due  to  the  fact  that  the  momentum  equation  is  not  at  steady  state . 

V.  LUMPED  PARAMETER  MODEL 

In  a  previous  paper, ^  a  lumped  parameter  model  was  derived  for  the 
injection  process.  The  density  is  assumed  to  be  constant  with  respect 
to  space  in  the  reservoir.  The  primary  assumption  is  that  the  mass  flux 
uA  is  also  constant  with  respect  to  space.  The  gradient  is  in  fact  very 
small,  since  the  liquid  is  not  very  compressible.  With  this  assumption, 
it  is  possible  to  integrate  the  momentum  equation  and  obtain 
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at 

The  momentum  equation  is  rearranged  so  the  area  terms  are  only  on  the 
left  side  of  the  equation  before  the  space  integration  is  performed. 

The  integral  of  the  inverse  of  the  cross  sectional  area  is  evaluated 
numerically.  In  the  code,  equation  (13)  is  integrated  in  parallel  with 
the  one -dimensional  equations. 

Like  many  lumped  parameter  models,  the  above  is  not  completely  self 
consistent.  Since  the  density  is  constant  in  the  reservoir,  so  is  the 
pressure  p^.  The  pressure  drop  is  idealized  as  occurring  instantan¬ 
eously  at  the  outflux  boundary.  Similarly,  the  assumption  that  the  mass 
flux  is  constant  implies  that  once  u^  is  known,  the  velocity  anywhere  in 
the  reservoir  can  be  determined  by  continuity.  Again,  the  change  in  the 
mass  flux  is  idealized  as  occurring  instantaneously  at  the  boundary. 

The  above  model  then  in  a  crude  fashion  implements  the  effects  of 
inertia,  as  there  is  a  delay  before  a  pressure  gradient  or  velocity 
gradient  has  an  effect  on  the  mass  flux  out  of  the  reservoir. 
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0.0087  ms  (line) 
0.0241  ms  (dot). 


Now  consider  why  the  pressure  wave  in  the  code  inverts  at  the  right 
hand  boundary.  The  pressure  at  the  boundary  is  kept  constant,  so  as  the 
high  pressure  point  approaches  the  boundary,  the  pressure  gradient 
becomes  steeper.  The  fluid  velocity  becomes  larger.  This  relatively 
large  velocity  causes  a  rarefaction  wave  to  form  and  propagate  back  into 
the  liquid  region. 

The  proper  boundary  conditions  for  the  outflow  boundary  are  not 
obvious.  A  number  of  workers  have  tried  to  develop  non- reflecting 
boundary  conditions.  If  the  boundary  is  artificial,  that  is,  introduced 
solely  to  reduce  the  size  of  the  region  of  integration,  pressure  waves 
should  exit  through  the  outflow  boundary.  However,  for  the  problems  of 
interest  here  non-reflecting  boundary  conditions  are  not  appropriate. 

At  the  boundary  there  is  an  abrupt  area  change  as  the  fluid  enters  the 
larger  combustion  chamber.  The  pressure  should  drop  rapidly  to  match 
the  combustion  chamber  pressure,  and  the  behavior  of  the  code  at  this 
boundary  is  reasonable . 

Next  consider  a  tube  where  the  radius  changes  rapidly  from  1  cm  to 
1/3  cm  between  x  -  -2.0  cm  and  x  -  -1.9  cm.  As  before,  a  pressure  wave 
is  started  at  the  left  boundary.  Figure  4  shows  the  pressure  profile 
after  the  original  wave  has  traveled  1.5  cm.  The  interesting  point  is 
that  the  change  in  area  splits  the  pressure  wave.  Part  of  the  wave 
continues  forward.  As  the  area  decreases,  the  pressure  increases. 
However,  part  of  the  wave  is  reflected  and  moves  back  toward  the  left. 

At  later  times  the  pressure  profile  becomes  very  complex  as  the  waves 
are  repeatedly  split. 

To  see  the  effects  of  an  sudden  increase  in  area,  consider  a  tube 
where  the  radius  changes  from  1/3  cm  to  1  cm  between  x  -  -2.0  cm  and  x  - 
-1.9  cm.  Figure  5  shows  the  profile  corresponding  to  Figure  4.  As 
before,  part  of  the  wave  continues  to  the  right.  As  the  area  Increases, 
the  pressure  decreases.  The  reflected  part  of  the  wave  is  inverted. 

The  procedure  is  that  discussed  above  for  the  right  hand  boundary. 


Figure  5. 


Pressure  Wave  in  a  Diverging  Cylinder 
Area  increases  at  x  -  -2.0  cm. 

Time  -  0.0087  ms. 
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VII.  INJECTION  FROM  A  CYLINDER 


Again  consider  a  cylinder  3  cm  long  with  a  radius  of  1  cm.  The 
back  wall  is  accelerated  smoothly  to  500  cm/s  in  the  first  0.5  ms.  The 
outflux  pressure  is  kept  fixed  at  10  MPa.  The  integration  is  carried 
out  to  1.0  ms.  Figure  6  shows  the  pressure  at  the  back  wall  as  a 
function  of  time.  As  the  back  wall  accelerates,  it  raises  the  pressure. 
There  is  a  time  delay  until  this  high  pressure  region  relaxes.  After 
the  piston  reaches  a  steady  velocity,  the  pressure  oscillates  around  the 
steady  state  value  of  10  MPa.  The  time  from  a  maximum  to  a  minimum 
value  of  the  pressure  matches  very  closely  the  time  for  a  sound  wave  to 
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ipace  profiles  of  the  pressure  at  the  later 
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Figure  7.  Injection  for  a  Cylinder.  Pressure  Profiles. 

Time  -  0.9333  ms  (line).  Time  -  0.9467  ms  (dot). 
Time  -  0.9602  ms  (dash). 


VIII.  TEST  PROBLEM 

As  a  more  realistic  test  problem,  an  orifice  with  a  typical  piston 
shape  is  considered,  but  without  a  central  bolt.  The  orifice  is  3  cm 
long,  and  the  area  change  from  left  to  right  is  a  factor  of  9  (see 
Figure  8).  As  above,  the  back  wall  is  accelerated  smoothly  to  500  cm/s 
in  the  first  0.5  ms.  The  outflux  pressure  is  kept  fixed  at  10  MPa.  Th< 
integration  is  carried  out  to  1.0  ms.  This  relatively  simple  problem 
will  be  used  to  compare  the  results  of  the  various  codes  and  for  grid 
checking . 
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Table  1  shows  a  set  of  runs  on  the  CRAY-XMP  for  the  one -dimensional 
code.  The  pressure  at  the  back  wall,  p^,  at  1.0  ms  is  given  as  the  most 
sensitive  output  value.  This  pressure  primarily  determines  the 
discharge  coefficient.  The  problem  has  essentially  converged  for  400 
grid  points.  The  problem  has  not  reached  steady  state  at  1.0  ms,  and 
the  discharge  coefficient  oscillates  around  the  value  one. 

Table  1.  Grid  Check  Runs  for  ODIN  Code. 


N 

Pi 

(MPa) 

CD 

Run  Time 

min: sec 

11 

11.76 

.861 

r 

1 

21 

12.44 

.822 

1 

41 

11.82 

.895 

2 

81 

11.73 

.917 

6 

161 

11.49 

.987 

22 

241 

11.49 

.986 

47 

321 

11.46 

.997 

1:22 

401 

11.42 

1.009 

2:07 

601 

11.44 

1.003 

4:41 

801 

11.44 

1.000 

8:17 

1001 

11.45 

1.007 

12:53 

The  test  problem  was  also  solved  using  the  MAGIC  code.  For 
simplicity  laminar  flow  was  assumed.  The  number  of  radial  grid  points 
NR  was  chosen  first.  The  axial  spacing  is  picked  so  that  the  grid 
quadrilaterals  are  square  in  the  first  straight  section  of  the 
reservoir.  The  quadrilaterals  in  the  last  straight  section  are 
rectangles  with  an  aspect  ratio  of  1.5.  The  axial  grid  size  decrease 
smoothly  in  the  converging  section.  This  choice  of  aspect  ratios  has 
been  observed  to  be  stable,  even  for  very  coarse  grids.  If  large  aspect 
ratios  are  chosen,  the  solution  will  break  up.  Table  2  shows  the  grid 


check  results.  The  numerical  convergence  is  very  slow,  and  the  solution 
has  not  converged  for  the  largest  grid  possible.  Part  of  the  problem  is 
the  upwind  differencing.  This  generates  a  numerical  viscosity  term  that 
depends  on  the  grid  spacing. ^  The  physical  viscosity  for  room 
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temperature  propellant  is  small  (kinematic  viscosity  -  0.05  stokes), 
and  the  numerical  viscosity  becomes  larger  than  the  physical  viscosity, 
even  for  a  fine  grid  for  this  simple  problem.  A  more  accurate  finite 
difference  method  appears  to  be  required  to  predict  the  discharge 
coefficients  for  the  more  complicated  problem  of  interest  (see  below) . 


Table  2. 

Grid  Check 

Runs  for 

MAGIC  Code. 

NR 

NX 

Pi 

CD 

Run  Time 

(MPa) 

hr  :min: sec 

6 

21 

14.94 

.540 

45 

11 

40 

13.47 

.645 

3:32 

21 

79 

12.73 

.726 

22:33 

41 

158 

12.32 

.788 

2:42:28 

61 

236 

12.16 

.816 

8:48:08 

81 

315 

12.07 

.834 

20:55:17 

The  MAGIC  code  is  partially  implicit,  and  should  not  be  restricted 
by  the  Courant  condition,  equation  (10).  However,  as  the  time  step  was 
increased,  the  code  did  many  iterations  on  the  implicit  pressure 
iteration  to  obtain  convergence ,  and  the  run  time  did  not  improve .  The 
code  was  finally  run  using  a  time  step  only  twice  the  Courant  limit. 

This  time  step  restraint  may  be  due  to  accuracy  considerations.  Since 
the  implicit  iterations  did  not  seem  to  be  very  helpful  for  this  type  of 
problem,  it  was  not  implemented  in  the  one -dimensional  code. 


Figure  8  shows  the  vector  plot  generated  by  the  MAGIC  code.  The 
solution  is  smooth  and  stable.  For  runs  using  realistic  piston  shapes, 
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such  phenomena  as  flow  separation  and  recirculation  have  not  been 
observed. 

It  is  more  interesting  to  plot  the  boundary  conditions  as  a 
function  of  time.  Figure  9  shows  the  influx  velocity  (specified  by  the 
user)  and  the  outflux  velocities  obtained  from  the  different  codes.  All 
three  codes  predict  outflux  velocities  very  close  to  the  steady  state 
values  (outflux  velocity  -  9  x  influx  velocity) .  A  blow  up  of  the 
outflux  velocities  (Figure  10)  shows  oscillations  due  to  the  pressure 
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Figure  9 . 


MAGIC  code  (line) 
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Figure  10. 


MAGIC  Code  (line).  ODIN  Code  (dot).  LPIN  (dash) 


Figure  11  shows  the  pressures  at  the  back  wall.  The  two-dimensional 
and  one -dimensional  codes  predict  qualitatively  very  similar  results. 

As  the  back  wall  accelerates,  it  raises  the  pressure.  There  is  a 
definite  time  delay  until  this  high  pressure  region  relaxes.  The 
differences  are  at  least  partly  due  to  the  lack  of  convergence  of  the 
MAGIC  code.  The  lumped  parameter  code  does  not  predict  the  initial 
large  overshoot  in  pressures .  Because  of  the  lack  of  spatial  resolution 
in  the  mass  flux,  LPIN  predicts  a  slightly  higher  mass  flux  out  of  the 
reservoir.  Even  a  slight  change  in  the  density  of  the  liquid  in  the 
reservoir  leads  to  a  large  pressure  difference.  There  is  a  secondary- 
effect  due  to  the  fact  that  LPIN  assumes  a  uniform  pressure  in  the 
reservoir  instead  of  a  pressure  gradient.  The  pressure  values  primarily 
determine  the  computed  discharge  coefficients  (Figure  12). 


Figure  11. 


time  (me) 


MAGIC  Code  (line).  ODIN  Code  (dot) 
LPIN  (dash)  . 


The  lumped  parameter  code  predicts  pressure  oscillations,  despite 
the  fact  that  there  are  no  pressure  waves  or  gradients  in  the  model. 
This  is  a  global  phenomenon.  Suppose  the  velocity  of  the  piston  is 
constant,  and  the  liquid  pressure  is  above  the  steady  state  pressure. 
The  fluid  will  be  accelerated.  As  the  fluid  is  injected  more  rapidly 
from  the  reservoir,  the  pressure  will  drop.  Since  the  model  involves  a 
transient  equation,  the  fluid  takes  time  to  slow  down,  and  will  be 
injected  more  rapidly  even  after  the  liquid  pressure  reaches  the  steady 
state  pressure.  So  the  liquid  pressure  will  end  up  below  the  steady 
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Figure  12. 
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MAGIC  Code  (line).  ODIN  Code  (dot). 
LPIN  (dash). 


The  one -dimensional  model  predicts  larger  oscillations.  These  no 
longer  correspond  to  the  time  scale  of  a  pressure  wave.  The  change  in 
area  causes  the  pressure  surges  to  be  broken  up  in  unpredictable  ways. 
These  pressure  surges  can  also  partially  cancel  one  another,  and  the 
damping  is  much  more  rapid  than  for  the  case  of  a  cylinder  without  area 
changes . 

The  MAGIC  code  predicts  results  very  similar  to  the  ODIN  code,  but 
the  oscillations  are  damped  more  rapidly.  This  is  at  least  partly  due 
to  the  large  numerical  viscosity.  It  may  also  be  partially  due  to  two- 
dimensional  effects  (such  as  wall  friction). 
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Data  from  the  BRL  30-mm  gun  has  been  analyzed. Here  the  data 
from  Round  8  is  selected.  The  chamber  pressure  and  the  piston  travel 
have  been  measured.  The  experimental  values  are  fitted  by  splines.  The 
chamber  pressure  fit  and  the  derivative  of  the  piston  travel  fit  are 
used  as  input  to  ODIN.  The  code  computes  the  back  wall  pressure  and 
hence  the  discharge  coefficients . 

The  liquid  pressure  at  the  back  wall  has  been  measured  using  mini- 
transducers .  There  has  been  some  concern  over  the  accuracy  of  these 
transducers.  A  partial  independent  check  is  possible.  The  injection  is 
caused  by  the  differential  piston.  At  steady  state,  the  liquid  pressure 
should  be  equal  to  the  chamber  pressure  times  the  hydraulic  difference 
(ratio  of  chamber  piston  area  over  liquid  reservoir  piston  area) . 

Figure  13  shows  a  comparison  between  the  measured  liquid  pressure  and 
the  chamber  pressure  times  the  hydraulic  difference.  The  two  curves  are 
in  reasonable  agreement  over  most  of  the  firing  cycle. 

The  transducer  block  (at  the  back  of  the  liquid  reservoir)  is 
mounted  on  Belleville  springs.  As  the  piston  begins  to  move,  it 
compresses  the  liquid,  which  pushes  back  on  the  transducer  block  and 
compresses  the  Belleville  springs.  This  helps  the  piston  to  clear  the 
O-ring  originally  sealing  the  reservoir  and  begin  the  injection  process. 
Very  little  injection  takes  place  until  the  springs  have  been 
compressed.  Rather  than  include  a  model  of  the  springs  in  the  inverse 
code,  the  assumption  is  made  that  the  piston,  liquid,  and  springs  move 
as  a  unit  until  the  springs  are  fully  compressed.  This  is  taken  as  time 
zero.  At  this  point  the  piston  is  assumed  to  begin  injecting  the 
liquid.  The  initial  conditions  are  not  quite  accurate,  since  the 
compression  of  the  springs  and  the  injection  process  are  not  actually 
separate  phases . 


tune  (ms) 


Figure  13. 


Figure  14  shows  the  central  bolt  and  piston  at  time  zero.  The  vent 
area  has  just  begun  to  open  up.  Figure  15  shows  the  outflux  velocity. 
The  models  predict  large  oscillations  near  the  start,  probably  due  to 
the  fact  that  the  initial  conditions  are  not  accurate.  After  this,  the 
predicted  one -dimensional  velocity  shows  very  good  agreement,  and  the 
predicted  zero -dimensional  velocity  is  still  reasonably  accurate. 
However,  the  agreement  for  the  pressures  is  not  as  good  (Figure  16). 
Between  1.0  and  2.5  ms,  both  models  predict  a  liquid  pressure  well  below 
the  values  measured.  This  leads  to  higher  discharge  coefficients 
(Figure  17).  The  one -dimensional  model,  like  the  lumped  parameter 


Figure  15. 
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Round  8  (line).  ODIN  Code 
(dot).  LPIN  Code  (dash). 


model,  shows  the  discharge  coefficient  increasing  rapidly,  rather  than 
the  slow  increase  derived  from  the  experimental  data.  The  one¬ 
dimensional  model  predicts  small  oscillations  in  the  liquid  pressure 
which  lead  to  large  oscillations  in  the  discharge  coefficient.  These 
oscillations  are  not  numerical  (they  are  independent  of  the  space  and 
time  discretization)  and  are  due  to  the  pressure  surges  being  split  by 
the  area  changes.  In  the  actual  flow  turbulence  and  friction  against 
the  walls  will  tend  to  damp  these  oscillations.  ODIN  and  LPIN  show  much 
better  agreement  for  the  liquid  pressure  here  than  for  the  simpler  test 
problem.  Both  the  back  wall  velocity  and  the  outflux  pressure  are  more 
complicated.  LPIN  alternately  overpredicts  and  underpredicts  the 
outflux.  The  errors  tend  to  cancel. 
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Figure  16.  Back  Wall  Pressures 
Round  8  (line). 

ODIN  Code  (dot). 
LPIN  Code  (dash). 


The  discharge  coefficient  is  a  useful  correlation  because  it  is  a 
measure  of  the  distance  from  steady  state  conditions.  The  pressure  and 
velocity  oscillations  at  the  start  of  the  firing  cycle  lead  to  very 
large  oscillations  in  the  discharge  coefficient.  A  system  without 
Belleville  springs  would  generate  cleaner  profiles,  since  the  Initial 
oscillations  in  pressure  and  velocity  would  not  be  generated. 
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Figure  17. 


2.0  2.5 

tine  (ms) 


Round  8  (line). 
ODIN  Code  (dot). 
LPIN  Code  (dash). 
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The  one -dimensional  code  does  not  include  the  effect  of  friction 

against  the  side  walls.  Friction  should  cause  a  higher  liquid  pressure 

and  damp  out  the  oscillations.  Since  this  could  lead  to  better 

agreement  with  the  experiment,  the  effect  of  friction  was  examined  using 

an  empirical  correlation.  For  fully  developed  turbulent  flow  in  a 

cylinder  of  diameter  D  and  length  L,  the  head  loss  due  to  frictional 
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effects  is  approximately 
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where  u  is  the  fluid  velocity  and  Re  is  the  Reynolds  number.  This  gives 
a  correction  to  the  steady  state  Bernoulli  equation.  The  head  loss 
leads  to  a  larger  pressure  difference  from  one  end  of  the  cylinder  to 
the  other.  The  correlation  is  applied  to  the  velocity  control  volume. 
For  the  diameter  of  the  cylinder  the  hydraulic  diameter  -  2.0  x  Gap 
is  substituted.  Adjusting  the  dimensions,  the  quantity 
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is  added  to  equation  (6) ,  the  finite  difference  form  of  the  momentum 
equation.  This  gives  an  estimate  to  the  effects  of  friction  across  the 
control  volume . 
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The  value  of  kinematic  viscosity  is  0.05  stokes.  Using  the  above 
correlation,  the  results  are  almost  identical  to  the  previous  answers. 
Even  if  the  kinematic  viscosity  is  increased  to  5.0  stokes,  the 
differences  are  minor  at  the  early  times  and  only  become  important  late 
in  the  cycle  (see  Figures  18-19).  While  the  above  approximation  is 
crude,  it  should  give  the  correct  order  of  magnitude  for  frictional 
effects.  The  tentative  conclusion  is  that  for  these  high  velocities  and 
relatively  low  viscosities,  frictional  effects  are  not  very  important. 
When  the  MAGIC  code  is  fully  operational,  this  conclusion  will  be 
checked. 
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Figure  19. 


X.  CONCLUSIONS 


Three  models  (lumped  parameter,  one -dimensional ,  and  two- 
dimensional)  have  been  developed  for  the  injection  of  liquid  propellant 
s  in  a  regenerative  gun.  For  a  simple  test  case,  the  models  show 

j  reasonable  agreement.  For  the  actual  gun  simulation,  only  the  two 

i  simpler  models  can  be  run.  They  both  predict  a  rapid  rise  to  roughly  a 

steady  state  discharge  coefficient.  The  results  do  not  agree  with 
values  for  the  discharge  coefficient  obtained  from  the  experimental 
data.  Preliminary  indications  are  that  the  two-dimensional  model  will 
give  similar  results.  However,  the  development  of  the  two-dimensional 
model  should  help  to  delineate  the  limitations  of  the  one -dimensional 
analysis . 

In  conclusion,  there  is  some  phenomenon  in  the  injection  process 
which  is  not  properly  taken  into  account.  The  modeling  of  the  liquid 
behavior  may  be  perfectly  adequate,  and  the  problems  due  to  other  parts 
of  the  system.  For  example,  the  chamber  pressure  is  assumed  to  be 
uniform,  and  the  pressure  measured  at  the  wall  is  taken  to  be  equal  to 
the  pressure  at  the  exit  of  the  reservoir.  This  is  a  boundary  condition 
for  the  injection  models  considered  here.  If  there  is  a  noticeable 
pressure  gradient  in  the  chamber,  this  condition  will  be  incorrect. 
Another  possibility  is  the  piston  motion.  The  velocity  of  the  piston  is 
another  boundary  condition,  and  the  position  of  the  piston  determines 
the  shape  of  the  liquid  reservoir.  If  the  piston  wobbles  about  the 
central  bolt,  this  will  effect  the  injection  area.  The  Belleville 
springs  may  have  a  more  complicated  behavior  than  than  assumed,  which 
would  effect  the  measured  position  of  the  piston.  Also,  the  piston 
travel  measurements  have  ignored  the  effect  of  the  recoil  of  the  gun. 
These  and  other  possible  contributions  to  the  piston  motion  are  being 
investigated. 
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LIST  OF  SYMBOLS 


speed  of  sound,  cm/s. 

o 

cross-sectional  area  of  liquid  reservoir,  cm  . 

Courant  number, 
discharge  coefficient, 
hydraulic  diameter,  cm. 

2 

conversion  constant,  10  g/s  -cm-MPa. 
bulk  modulus ,  MPa . 

bulk  modulus  of  the  liquid  at  zero  pressure,  MPa. 

derivative  of  the  bulk  modulus. 

mass  in  a  region,  g. 
pressure,  MPa. 

Reynolds  number  (u  D^/v) . 
surface  area,  cm^ . 
time ,  s . 

fluid  velocity,  cm/s. 


grid  velocity,  cm/s. 


I  li 


3 

V  volume ,  cm  . 

v  kinematic  viscosity,  stokes. 

3 

pQ  density  of  the  liquid  at  atmospheric  pressure,  g/cm  . 

3 

P  density  of  the  liquid  in  the  reservoir,  g/cm  . 
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